//
//  Jacobian.cpp
//  EMT_C
//
//  Created by Duong Ngo on 8/27/16.
//  Copyright © 2016 Duong Ngo. All rights reserved.
//

#include "Jacobian.hpp"

#include <math.h>

void find_jacobian_structure(const std::map<std::string,int>& f, Index* iRow, Index* jCol, int& dem)
{
    //-------------------------------------------------------------
    // BANKERS
    //-------------------------------------------------------------
    //g0
    int pg=0;
    iRow[dem]=pg; jCol[dem]=f.at("s.gamma"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.cb"); ++dem;
    
    
    //g1
    ++pg;
    iRow[dem]=pg; jCol[dem]=f.at("s.Rf"); ++dem;
    
    //g2
    ++pg;
    iRow[dem]=pg; jCol[dem]=f.at("s.gamma"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.Rm"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.mur"); ++dem;
    
    //g3
    ++pg;
    iRow[dem]=pg; jCol[dem]=f.at("s.gamma"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.mur"); ++dem;

    //g4
    ++pg;
    iRow[dem]=pg; jCol[dem]=f.at("s.ql"); ++dem;

    //g5
    ++pg;
    iRow[dem]=pg; jCol[dem]=f.at("s.tau"); ++dem;

    //g6
    ++pg;
    iRow[dem]=pg; jCol[dem]=f.at("s.m"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.Rm"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.ql"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.s"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.bh"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.cb"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.n"); ++dem;

    //g7
    ++pg;
    iRow[dem]=pg; jCol[dem]=f.at("s.n"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.m"); ++dem;
    
    //g8
    ++pg;
    iRow[dem]=pg; jCol[dem]=f.at("s.n"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.bh"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.m"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.muc"); ++dem;
    
    
    //g9
    ++pg;
    iRow[dem]=pg; jCol[dem]=f.at("s.bh"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.s"); ++dem;

   
    //-------------------------------------------------------------
    // HOUSEHOLDS
    //-------------------------------------------------------------
    
    //g10
    ++pg;
    iRow[dem]=pg; jCol[dem]=f.at("s.lama"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.etaz"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.ch"); ++dem;
 
    //g11
    ++pg;
    iRow[dem]=pg; jCol[dem]=f.at("s.lamb"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.ch"); ++dem;
    
    //g12
    ++pg;
    iRow[dem]=pg; jCol[dem]=f.at("s.lama"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.Rm"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.lamb"); ++dem;

    //g13
    ++pg;
    iRow[dem]=pg; jCol[dem]=f.at("s.ql"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.lamb"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.etab"); ++dem;

    //g14
    ++pg;
    iRow[dem]=pg; jCol[dem]=f.at("s.lamb"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.pm"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.lama"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.k"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.y"); ++dem;


    //g15
    ++pg;
    iRow[dem]=pg; jCol[dem]=f.at("s.Rm"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.m"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.ql"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.s"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.ch"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.i"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.bh"); ++dem;
    
    
    //g16
    ++pg;
    iRow[dem]=pg; jCol[dem]=f.at("s.bh"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.etab"); ++dem;
    
    //g17
    ++pg;
    iRow[dem]=pg; jCol[dem]=f.at("s.pm"); ++dem;

    //g18
    ++pg;
    iRow[dem]=pg; jCol[dem]=f.at("s.y"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.k"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.l"); ++dem;

    
    //g19
    ++pg;
    iRow[dem]=pg; jCol[dem]=f.at("s.u"); ++dem;

  
    //g20
    ++pg;
    iRow[dem]=pg; jCol[dem]=f.at("s.y"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.cb"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.ch"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.i"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.bh"); ++dem;
    
    //g21
    ++pg;
    iRow[dem]=pg; jCol[dem]=f.at("s.k"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.i"); ++dem;
    
    //Inflation
    
    //g22
    ++pg;
    iRow[dem]=pg; jCol[dem]=f.at("s.ip"); ++dem;

    //g23
    ++pg;
    iRow[dem]=pg; jCol[dem]=f.at("s.pie"); ++dem;
    
    
    //g24
    ++pg;
    iRow[dem]=pg; jCol[dem]=f.at("s.l"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.pm"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.y"); ++dem;
    iRow[dem]=pg; jCol[dem]=f.at("s.lama"); ++dem;

    
    
}
//***********************************************************************************

//***********************************************************************************

void find_jacobian_values(const variables& s, Number* values, int& dem)
{
    //-------------------------------------------------------------
    // BANKERS
    //-------------------------------------------------------------
    //g0
    int pg=0;
//    iRow[dem]=pg; jCol[dem]=f.at("s.gamma"); ++dem;
    values[dem]= s.cb;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.cb"); ++dem;
    values[dem]= s.gamma;
    ++dem;
    
    
    //g1
    ++pg;
//    iRow[dem]=pg; jCol[dem]=f.at("s.Rf"); ++dem;
    values[dem]= 1;
    ++dem;
    
    //g2
    ++pg;
//    iRow[dem]=pg; jCol[dem]=f.at("s.gamma"); ++dem;
    values[dem]= 1- beta_b*s.Rm;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.Rm"); ++dem;
    values[dem]= -beta_b*s.gamma;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.mur"); ++dem;
    values[dem]= -varphi;
    ++dem;
    
    //g3
    ++pg;
//    iRow[dem]=pg; jCol[dem]=f.at("s.gamma"); ++dem;
    values[dem]= 1-beta_b*Rn;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.mur"); ++dem;
    values[dem]= -1;
    ++dem;
    
    //g4
    ++pg;
//    iRow[dem]=pg; jCol[dem]=f.at("s.ql"); ++dem;
    values[dem]= 1;
    ++dem;
    
    //g5
    ++pg;
//    iRow[dem]=pg; jCol[dem]=f.at("s.tau"); ++dem;
    values[dem]= 1;
    ++dem;
    
    //g6
    ++pg;
//    iRow[dem]=pg; jCol[dem]=f.at("s.m"); ++dem;
    values[dem]= 1-s.Rm;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.Rm"); ++dem;
    values[dem]= -s.m;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.ql"); ++dem;
    values[dem]= -s.s;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.s"); ++dem;
    values[dem]= -s.ql;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.bh"); ++dem;
    values[dem]= deltab- thetaa;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.cb"); ++dem;
    values[dem]= -1;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.n"); ++dem;
    values[dem]= (Rn-1);
    ++dem;
    
    //g7
    ++pg;
//    iRow[dem]=pg; jCol[dem]=f.at("s.n"); ++dem;
    values[dem]= 1;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.m"); ++dem;
    values[dem]= -varphi;
    ++dem;
    
    //g8
    ++pg;
//    iRow[dem]=pg; jCol[dem]=f.at("s.n"); ++dem;
    values[dem]= 1;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.bh"); ++dem;
    values[dem]= 1-kappa;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.m"); ++dem;
    values[dem]= -1;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.muc"); ++dem;
    values[dem]= -1;
    ++dem;
    
    
    //g9
    ++pg;
//    iRow[dem]=pg; jCol[dem]=f.at("s.bh"); ++dem;
    values[dem]= (1-deltab);
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.s"); ++dem;
    values[dem]= -1;
    ++dem;
    
    
    //-------------------------------------------------------------
    // HOUSEHOLDS
    //-------------------------------------------------------------
    
    //g10
    ++pg;
//    iRow[dem]=pg; jCol[dem]=f.at("s.lama"); ++dem;
    values[dem]= -1;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.etaz"); ++dem;
    values[dem]= -1;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.ch"); ++dem;
    values[dem]= -1/s.ch/s.ch;
    ++dem;
    
    //g11
    ++pg;
//    iRow[dem]=pg; jCol[dem]=f.at("s.lamb"); ++dem;
    values[dem]= s.ch;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.ch"); ++dem;
    values[dem]= s.lamb;
    ++dem;
    
    //g12
    ++pg;
//    iRow[dem]=pg; jCol[dem]=f.at("s.lama"); ++dem;
    values[dem]= 1;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.Rm"); ++dem;
    values[dem]= -beta_h*s.lamb;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.lamb"); ++dem;
    values[dem]= -beta_h*s.Rm;
    ++dem;
    
    //g13
    ++pg;
//    iRow[dem]=pg; jCol[dem]=f.at("s.ql"); ++dem;
    values[dem]= s.lamb- beta_h*deltab*s.lamb;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.lamb"); ++dem;
    values[dem]= s.ql- beta_h*(deltab+ deltab*s.ql);
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.etab"); ++dem;
    values[dem]= -pe*2*s.etab;
    ++dem;
    
    //g14
    ++pg;
//    iRow[dem]=pg; jCol[dem]=f.at("s.lamb"); ++dem;
    values[dem]= 1- beta_h*(1-deltak);
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.pm"); ++dem;
    values[dem]= -beta_h*alphaa*s.lama*s.y/s.k;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.lama"); ++dem;
    values[dem]= -beta_h*alphaa*s.pm*s.y/s.k;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.k"); ++dem;
    values[dem]= -beta_h*alphaa*s.pm*s.lama*s.y*(-1)/s.k/s.k;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.y"); ++dem;
    values[dem]= -beta_h*alphaa*s.pm*s.lama/s.k;
    ++dem;


    
    //g15
    ++pg;
//    iRow[dem]=pg; jCol[dem]=f.at("s.Rm"); ++dem;
    values[dem]= s.m;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.m"); ++dem;
    values[dem]= s.Rm;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.ql"); ++dem;
    values[dem]= s.s;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.s"); ++dem;
    values[dem]= s.ql;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.ch"); ++dem;
    values[dem]= -1;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.i"); ++dem;
    values[dem]= -1;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.bh"); ++dem;
    values[dem]= -deltab;
    ++dem;
    
    
    //g16
    ++pg;
//    iRow[dem]=pg; jCol[dem]=f.at("s.bh"); ++dem;
    values[dem]= -1;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.etab"); ++dem;
    values[dem]= 1;
    ++dem;
    
    //g17
    ++pg;
//    iRow[dem]=pg; jCol[dem]=f.at("s.pm"); ++dem;
    values[dem]= 1;
    ++dem;
    
    //g18
    ++pg;
//    iRow[dem]=pg; jCol[dem]=f.at("s.y"); ++dem;
    values[dem]= 1;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.k"); ++dem;
    values[dem]= -alphaa*pow(s.k,alphaa-1)*pow(s.l,1-alphaa);
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.l"); ++dem;
    values[dem]= -(1-alphaa)*pow(s.l,-alphaa)*pow(s.k,alphaa);
    ++dem;
    
    //g19
    ++pg;
//    iRow[dem]=pg; jCol[dem]=f.at("s.u"); ++dem;
    values[dem]= 1;
    ++dem;
    
    
    //g20
    ++pg;
//    iRow[dem]=pg; jCol[dem]=f.at("s.y"); ++dem;
    values[dem]= 1;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.cb"); ++dem;
    values[dem]= -1;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.ch"); ++dem;
    values[dem]= -1;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.i"); ++dem;
    values[dem]= -1;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.bh"); ++dem;
    values[dem]= -thetaa;
    ++dem;
    
    //g21
    ++pg;
//    iRow[dem]=pg; jCol[dem]=f.at("s.k"); ++dem;
    values[dem]= deltak;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.i"); ++dem;
    values[dem]= -1;
    ++dem;
    
    //Inflation
    
    //g22
    ++pg;
//    iRow[dem]=pg; jCol[dem]=f.at("s.ip"); ++dem;
    values[dem]= 1;
    ++dem;
    
    //g23
    ++pg;
//    iRow[dem]=pg; jCol[dem]=f.at("s.pie"); ++dem;
    values[dem]= 1;
    ++dem;
    
    
    //g24
    ++pg;
//    iRow[dem]=pg; jCol[dem]=f.at("s.l"); ++dem;
    values[dem]= chi*(nu+1)*pow(s.l, nu);
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.pm"); ++dem;
    values[dem]= -(1-alphaa)*s.y*s.lama;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.y"); ++dem;
    values[dem]= -(1-alphaa)*s.pm*s.lama;
    ++dem;
//    iRow[dem]=pg; jCol[dem]=f.at("s.lama"); ++dem;
    values[dem]= -(1-alphaa)*s.pm*s.y;
    ++dem;
}























